Numerical Improvement of the Discrete Element 
Method applied to Shear of Granular Media 
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We present a detailed analysis of the bounds on the integration step in Discrete Ele- 
ment Method (DEM) for simulating collisions and shearing of granular assemblies. 
We show that, in the numerical scheme, the upper limit for the integration step, usu- 
ally taken from the average time tc of one contact, is in fact not sufficiently small to 
guarantee numerical convergence of the system during relaxation. In particular, we 
study in detail how the kinetic energy decays during the relaxation stage and com- 
pute the correct upper limits for the integration step, which are significantly smaller 
than the ones commonly used. In addition, we introduce an alternative approach, 
based on simple relations to compute the frictional forces, that converges even for 
integration steps above the upper limit. 



1 Introduction 

One of the standard approaches to model the dynamics of granular media is to use 
the Discrete Element Method (DEM) [1,2, 3], e.g. to study shear [4, 5, 6, 7, 8]. Some 
problems may arise due to the need to use large integration steps to perform numer- 
ical simulations with reasonable computational effort, without compromising the 
overall convergence of the numerical scheme. For slow shearing, the convergence 
of the numerical schemes is particularly important when studying for instance the 
occurrence of avalanches [8] and the emergence of ratcheting in cyclic loading [9]. 
Usually, one assumes an upper limit for the admissible integration steps, based on 
empirical reasoning [10]. 

The aim of this paper is twofold. First, we show that an integration step able to 
guarantee convergence of the numerical scheme, must in general be smaller than a 
specific upper limit, significantly below the commonly accepted value [5, 10, 11]. 
This upper limit strongly depends on (i) the accuracy of the approach used to cal- 
culate frictional forces between particles, (ii) on the corresponding duration of the 
contact and (iii) on the number of degrees of freedom. Second, we address the spe- 
cific case of slow shearing, for which the above limit is too small to allow for rea- 
sonable computation time. To overcome this shortcoming we propose an alternative 
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Fig. 1. Illustration of two overlapping particles. The overlap region A between particles fully 
characterizes the contact force F"^ (see text). 

approach that corrects the frictional contact forces, when large integration steps are 
taken. In this way, we enable the use of considerably larger integration steps, assur- 
ing at the same time the convergence of the integration scheme. 

We start in Sec. 2 by presenting in some detail the Discrete Element Method [1, 
10, 12]. Sections 3 and 4 describe respectively the dependence on the integration 
step and the improved algorithm. Discussions and conclusions are given in Sec. 5. 



We consider a two-dimensional system of polygonal particles, each one having two 
linear and one rotational degree of freedom. The evolution of the system is given 
by Newton's equations of motion, where the resulting forces and moments acting 
on each particle i are given by the sum of all forces and momenta applied on that 
particle: 



where rm denotes the mass of particle i, U its moment of inertia, and I the branch 
vector which connects the center of mass of the particle to the application point of 
the contact force or boundary force F\. The sum in c is over all the particles 
in contact with polygon i, and the sum in ct is over all the vertices of polygon i in 
contact with the boundary. One integrates Eqs. (1) for all particles i = 1, . . . ,N and 
obtains the evolution of their centers of mass and rotation angles 9i. 



2 The model 
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Further, during loading, particles tend to deform each other This deformation of 
the particles is usually reproduced by letting them overlap [1, 12], as illustrated in 
Fig. 1. The overlap between each pair of particles is considered to fully characterize 
the contact. Namely, the normal contact force is assumed to be proportional to the 
overlap area [17] and its direction perpendicular to the plane of contact, which is 
defined by the intersection between the boundaries of the two particles. 

All the dynamics is deduced from the contact forces acting on the particles. The 
contact forces, F'^, either between particles or with the boundary, are decomposed 
into their elastic and viscous contributions, F'^ and F"" respectively, yielding F^ = 

The viscous force is important for maintaining the numerical stability of the 
method and to take into account dissipation at the contact. This force is calculated 
as [12] 

= -rurvv'', (2) 

where m,. = (l/m^ + l/mj)^^ is the reduced mass of the two particles, i and j, 
and V is the damping coefficient. 

The elastic part of the contact force is what will be carefully studied, since it 
is what determines the accuracy of the integration scheme. The term F''' is simply 
given by the sum of the normal and the tangential components, with respect to the 
contact plane, namely 

F" ^ F^^n'' + F^i\ (3) 
where the normal component reads 

F'^ = -KA/l,, (4) 

with fc„ the normal stiffness, A the overlap area and Ic ^ ri + rj the characteristic 
length of the contact between particles i and j, with = Aij^-n and Ai the area 
of the particle i (and similarly for particle j). 

Using an extension of the Cundall-Stack spring [12], which considers the tan- 
gential force to be proportional to the elastic elongation ^ of an imaginary spring at 
the contact, one defines the static frictional force between each pair of particles in 
contact, as 

Ft = -kti, (5) 

where kt is the tangential stiffness. This tangential force assumes each contact as 
being described by a damped oscillator with some frequency uj (see below). 
The elastic elongation ^ in Eq. (5) is updated as 

e(t + At) = + VtAt (6) 

where At is the time step of the DEM simulation, and v'j: is the tangential component 
of the relative velocity 



v'^ = Vi — Vj + Wi X li — Wj X Ij , 



(7) 
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at the contact point. Here, Vj and Vj are the linear velocities of the centers of mass 
and Wi and Wj the angular velocities of the particles around the corresponding 
centers of mass. 

The tangential elastic elongation ^ changes according to Eq. (6) whenever the 
condition \F*\ < fj,F" is satisfied, whereas, when the Coulomb limit condition 

I = fiF" is reached, sliding is enforced by keeping constant the tangential force 
Fg and assigning to ^ its extreme values ±jiknA/ (ktlc)- This latter Coulomb condi- 
tion corresponds to the regime where particles behave inelastically, while the former 
inequality describes the forces when the particles behave elastically. Parameter /^t is 
the inter-particle friction coefficient. 

In DEM, one of the numerical integration schemes usually used to solve the 
equations of motion above is Gear's predictor-corrector scheme [10]. This scheme 
consist of three main stages, namely prediction, evaluation and correction. 

In the prediction stage the linear and angular positions, velocities and higher- 
order time derivatives are updated by expansions in Taylor series using the current 
values of these quantities [10, 13]. In particular, one extracts a predicted position 
rP (t + At) and acceleration (t + At) for the center of mass of a given particle and 
the predicted angular displacement 9^ {t + At) and angular acceleration 6 {t + At) 
of that particle around its center of mass. 

During the evaluation stage, one uses the predicted coordinates to determine 
the contact force F1_^_^f at time t + At. Since the method is not exact, there is a 
difference between the acceleration r{t + At) — F'j:_^^^/m and the value obtained 
in the prediction stage, namely 

Ar = r{t + At)-rP{t + At). (8) 

The difference Ar in Eq. (8) is then finally used in the corrector step to cor- 
rect the predicted position and its time derivatives, using proper weights for each 
time derivative [10], that depend upon the order of the algorithm and the differential 
equation being solved. These corrected values are the ones used for the next inte- 
gration step t + At. This same procedure is also applied to the rotation angles 6i 
around the center of mass as well as to their time derivatives, yielding the correction 
AO. 

In our simulations we integrate equations of the form r = f{r, f), using a fifth 
order predictor-corrector algorithm that has a numerical error proportional to (At)^ 
for each integration step [10]. However, as will be seen in Sec. 3, (At)^ is not the 
numerical error of the full integration scheme, since Eq. (5), used to calculate the 
frictional force, is of order (At)^. 

For a certain value of normal contact stiffness fc„, almost any value for the nor- 
mal damping coefficient f„ might be selected. Their relation defines the restitution 
coefficient e obtained experimentally for various materials [14]. The restitution co- 
efficient is given by the ratio between the relative velocities after and before a colli- 
sion. In particular, the normal restitution coefficient e„ can be written as a function 
of fc„ and i^n [15], namely 




Fig. 2. (Color online) Sketch of the system of 256 particles (green particles) under shearing 
of top and bottom boundaries (blue particles). Horizontally periodic boundary conditions are 
considered and a constant low shear rate is chosen (see text). 



where cu = y^ojQ — if is the frequency of the damped oscillator, with lliq = 
yjkn/mr the frequency of the elastic oscillator, nir the reduced mass and ?/ = 
Vn/{2mr) the effective viscosity. The tangential component et of the restitution co- 
efficient is defined similarly using kt and vt in Eq. (9). 

Therefore, a suitable closed set of parameters for this model are the ratios kt / fc„ 
and et / e„ (or vt / f„), together with the normal stiffness fc„ and the interparticle 
friction /^t. 

The entire algorithm above relies on a proper choice of the integration step At, 
which should neither be too large to avoid divergence of the integration nor too small 
avoiding unreasonably long computational time. The determination of the optimal 
integration step varies from case to case and there are two main criteria to estimate 
an upper bound for admissible integration steps. 

The first criterion is to use the characteristic period of oscillation [10], defined 

as 
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t,=2wJ^, (10) 



where (to) is the smallest particle mass in the system. For a fifth order predictor- 
corrector integration scheme, it is usually accepted that a safe integration step should 
be below a threshold of At < ts/10 [10]. 

The second criterion is to extract the threshold from local contact events [5, 1 1, 
15], namely from the characteristic duration of a contact: 

TT 

tc^ , =■ (11) 

Typically tc — ts/2, and therefore in such cases, one considers an admissible inte- 
gration step as At < tc/5 [15, 16]. 

In the next section we will study in detail the integration for different values of 
the model parameters. 



3 The choice of the integration step 

We simulate the relative motion of two plates shearing against each other [7, 8, 17]. 
Considering a system of 256 particles as illustrated in Fig. 2, where both top and 
bottom boundaries move in opposite directions with a constant shear rate 7. The 
top and bottom layer of the sample have fixed boundary conditions, while horizon- 
tally we consider periodic boundary conditions. The volumetric strain is suppressed, 
i.e. the vertical position of the walls is fixed and there is no dilation. The particles of 
the fixed boundary are not allowed to rotate or move against each other The shear 
rate is 7 = 1.25 • 10~^s~^, the parameter values are kn ~ 400 N/m, e„ = 0.9875, 
and /i = 0.5. The relation kt/kn is chosen such that kt < kn, similarly to previ- 
ous studies [12, 15, 18], namely kt/kn — 1/3. Further, for simplicity we consider 
i^t/i^n — kt/kn, which when substituted in Eq. (9) yields et/e„ = 1.0053. 

By integrating such a system of particles using the scheme described in the pre- 
vious section, one can easily compute the kinetic energy of a given particle i, 

Ek{^) = l{m,r',+I^iJ!), (12) 

where velocity r is computed from the predictor-corrector algorithm, li is the mo- 
ment of inertia of the polygon and uji is the angular velocity. 

In Fig. 3 we show the evolution of the kinetic energy for two different At = 
0.001 s and 0.005 s. As one sees, frictionless particles (Fig. 3a and 3c) have an Ek 
that does not change for different integration steps, while when /i = 0.5 (Fig. 3b and 
3d) the evolution of Ek strongly depends on At. In Fig. 3e we plot the cumulative 
difference AEk between the values of Ek taken for each integration step. Here, one 
sees that in the absence of friction AEk is significantly lower then when friction is 
present. 
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Fig. 3. Dependence of the numerical scheme on the integration step At and the friction co- 
efficient fj,, by plotting the kinetic energy Ek as a function of time, for (a) At — 10^'^ s and 
= (no friction), (b) At = 10"^ s and = 0.5, (c) At ^ 5 x lO"'^ s and = 0, and 
(d) Z\t = 5 X 10"=* s and ^ = 0.5. In (e) we show the difference between the values of Ek 
obtained with the two values of At. Here, k„ = 400 and the parametric relations in Eq. (10) 
and (11) are used (see text). 



The two values of At used in Fig. 3 can be written as At = 13/500tc and 
13/2500f:c. Thus, we conclude that the expected upper limit ~ tc/10 is still too large 
to guarantee convergence of the integration scheme when friction is considered. 
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Fig. 4. Sketch of the stress controlled test of two particles (discs). The particle located at Vi 
remains fixed, while the particle at rj is initially touching particle i. The vector Ri j connect- 
ing the center of mass of particles i and j is initially oriented 45° with respect to the x-axis. 
After applying the constant force F to disc j, the system relaxes to a new position (dashed 
circumference). Between its initial and final position particle j undergoes a displacement Ar 
and a rotation A9 (see text). 




Fig. 5. The relaxation of the system of two discs sketched in Fig. 4. Here we plot the kinetic 
energy Ek as a function of time t (in units of tc) for different integration increments At and 
using a stiffness k„ — 4 x 10* N/m and a friction coefficient fi — 500. The large fi value 
is chosen so that the system remains in the elastic regime. As one sees, the relaxation time 
tj{ converges to a constant value when At is sufficiently small (see text). This discrepancy 
between the values of tn when different integration steps are used does not occur in the 
absence of friction (jj, = 0), as illustrated in the inset. The slope of the straight lines is 
-1/tR (see Eq. (13)). 
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Next we will perform a careful analysis to obtain a proper integration step as 
function of the parameters of our model. For that, we consider the simple situation 
of two circular particles and study the kinetic energy of one of them under external 
forcing, as sketched in Fig. 4. We start with two touching discs, i and j, where one 
of them, say i, remains fixed, while the other is subject to a force F perpendicular 
to its surface (no external torque is induced) along the .x-axis. As a result of this 
external force, the disc j undergoes rotation. The contact force is obtained from the 
corresponding springs that are computed, as described in Sec. 2, and acts against the 
external force. This results in an oscillation of disc j till relaxation (dashed circle 
in Fig. 4) with a final center of mass displacement of AR and a rotation around the 
center of mass of A6. Since F is kept constant, the procedure is stress controlled. 

For the two discs sketched in Fig. 4, we plot in Fig. 5 the kinetic energy as 
a function of time, from the beginning until relaxation, using different integration 
steps, namely At ^ 10~\ 2 x 10"\ 10"^, 10"^ and lO"'^ s in units of tc. As we 
see, the kinetic energy decays exponentially, 

£;,(0=i?rexp(-^), (13) 

where tfl is a relaxation time whose value clearly depends on the integration step 
At. As illustrated in the inset of Fig. 5, this change in t^ is not observed when 
friction is absent (/i = 0), since no tangential forces are considered (F^ = 0). 
Next, we will show that this dependence of t r on At vanishes for 

At<Tt{kn,fl)tc, (14) 

where Tt{kn, /i) is a specific function that is determined below. Notice that, in our 
case, the only free parameters on which Tt may depend are the friction fi and the 
normal stiffness fc„, since we consider a fixed restitution coefficient in the normal 
direction, e„ — 0.9875 and fixed relations fcf/fc„ = 1/3 and et/e„ = 1.0053. 

Figure 6a shows the relaxation time tj^ of the kinetic energy of the two-particle 
system for different values of stiffnesses, namely for fc„ = 1,50,200,10^ and 
10^ N/m. For all fc„ values, one sees that, for decreasing At, the relaxation time 
t]^ increases until it converges to a maximum. The stabilization of t^ occurs when 
At is small compared to the natural period l/wg of the system. We define Tt as the 
largest value of At for which we have this maximal relaxation time. 

As shown in Fig. 6b, all curves in 6a can be collapsed by using the normalized 
integration step At/tc- From Eq. (1 1) we calculate the contact times corresponding 
to these fc„ values as t^ 1.969, 0.278, 0.139, 9.8 x 10"^ and 9.8 x 10"^ s respec- 
tively. In fact, as shown in the inset of Fig. 6b the relaxation time scales with the 
stiffness as tc ~ kn'^^^ (see Eq. (11)). 

From Fig. 6 one can conclude that the relaxation time converges when the in- 
tegration step obeys Eq. (14) with Tt = 10~^ (dashed vertical line in Fig. 6b). We 
simulate the system also for ^ ~ 0.005, 0.005, 0.05, 0.5, 5, 50 and 500 and similar 
results were obtained. 

In Fig. 6 both translation and rotation of particles are considered. The rotation 
of particles is usually of crucial interest, for instance to simulate rolling [7, 19]. 
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Fig. 6. The relaxation time tR (in units of tc) as a function of (a) the integration step At and 
(b) the normalized integration step At/tc, where the contact time tc is defined in Eq. (11). 
Here the friction coefficient is kept fixed /i — 500 and different stiffnesses fc„ (in units of 
N/m) are considered. The quotient At/tc collapses all the curves for different k„. We find 
tc ^ kn as illustrated in the inset (see Eq. (11)). As a final result one finds a constant 
Tt = 10""^ (dashed vertical line). For other values of the friction coefficient one observes 
similar results. 

But, suppressing rotation can also be of interest for instance, when simulating fault 
gouges: by hindering the rotation of particles, one can mimic young faults where a 
strong interlocking between the constituent rocks is expected [7]. 

To study this scenario, we present in Fig. 7a the relaxation time for the same 
parameter values as in Fig. 6, now disabling rotation. Here, we obtain a constant 
Tt = 10^^ also, independent of fc„, one order of magnitude smaller than the previ- 
ous value in Fig. 6. In other words, when rotation is suppressed, one must consider 
integration steps typically one order of magnitude smaller than in the case when the 
discs are able to rotate. This can be explained as follows. 

When suppressing rotation, one restricts the system to have a single degree of 
freedom. All energy stored in the rotational degree of freedom through the integra- 
tion of the equations of motion is suppressed. This effectively acts like an increase 
of the friction coefficient, making the system more sensitive to the integration step, 
i.e. yielding a smaller bound Tt- By comparing Fig. 7a with Fig. 6a, one sees that 
the relaxation time t b, is smaller when rotation is suppressed. 

From the bounds on the integration steps obtained above, one realizes that, in 
general, the correct integration step must be significantly smaller than usually as- 
sumed. 
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Fig. 7. The relaxation time tn (in units of tc) of the kinetic energy as a function of the 
normalized integration step At/tc, when rotation is suppressed, (a) n = 500 and different 
values of k„ and for (b) fc„ = 4 x 10* and different values of fi. The dashed horizontal line 
/i = in (b) indicates the relaxation time of the kinetic energy in the absence of friction (see 
text). 
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Fig. 8. The relaxation time tn (in units of tc) as a function of the friction coefficient /i 
when rotation is suppressed. Here fc„ = 4 x 10* N/m which corresponds to a contact time 
tc ~ 9.8 X 10^^ s. The normalized integration step is At/tc = 10^^. 
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While Fig. 7a clearly shows that tfj does not depend on the stiffness fc„, from 
Fig. 7b one sees that the same is not true for the friction coefficient /i. Indeed, from 
Fig. 8 one sees that there is a change of the relaxation time around ^ — 1. Here, the 
values correspond to a normalized integration step At/tc = 10~^ for which tn has 
already converged. This might be explained by considering the fact that for large 
values of ^ the contact is essentially non-sliding, which induces a faster relaxation 
than for smaller ^ values. 

It is important to stress that all the results above were taken within the elastic 
regime, since the dependence on At does not occur when the Coulomb condition 
is fulfilled (inelastic regime). This fact indicates that the improvements in the algo- 
rithm should be implemented when computing the elastic component of the tangen- 
tial contact force, in Eq. (6), as explained in the next Section. 

4 Improved approach to integrate the contact force 

In this Section we will describe a technique to overcome the need of very small inte- 
gration steps. As shown previously, when using Cundall's spring[4], the relaxation 
time of the two particles only converges when At is a small fraction Ti of the contact 
time tc- This is due to the fact that the elastic elongation is assumed to be linear in 
At, i.e. the finite difference scheme in Eq. (6) is very low order, {AtY, compromis- 
ing the convergence of the numerical scheme that is of order [AtY" . Therefore, the 
most plausible way to improve our algorithm is by choosing a different expression 
to compute the elastic tangential elongation ^ without using Eq. (6). 

We will introduce an expression for ^ that contains only the quantities computed 
in the predictor step. In this way we guarantee that ^ has errors of the order of {At)^, 
instead of [At)^, as it is the case of Eq. (6). Let us illustrate our approach on the 
simple system of two discs considered in the previous Section (see Fig. 4). 

On one side, if rotation is not allowed, the elastic elongation ^ depends only on 
the relative position of the two particles. In this case we substitute Eq. (6) by the 
expression 

+ At) = (0 + ^^{m,{t + At) m^{t)) ■ ds) 

where and aj are the radii of the discs i and j respectively, i?y is the vector 
joining both centers of mass and points in the direction i j (see Fig. 4). Index p 
indicates quantities derived from the coordinates computed at the predictor step. 

On the other side, if iiy is kept constant and only rotation is possible, particle 
j will have an elongation ^ that depends only on its rotation between time t and the 
predictor step: 

£![°'\t + At) = £}[°'\t) + iff^it + At) - e^{t))a,, (16) 

where Qp and 9{t) are the angles of some reference point on particle j at the predictor 
step and at time t respectively. 
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Fig. 9. The relaxation time (in units of tc) using Eqs. (15) and (16) between two discs, as 
illustrated in Fig. 4. For the three cases when considering only rotation, only translation or 
both, the relaxation time remains constant independent of the integration step. 

When both translation and rotation of particle j occur, then the elongation is the 
superposition of both contributions, yielding = + 

Figure 9 shows the relaxation time tn as a function of the integration step for the 
three situations above, namely when only rotation is considered, when only trans- 
lation is considered and when both rotation and translation are allowed. As we see 
for all these cases, the relaxation time is independent on the integration step. This is 
due to the fact that all quantities in the expression for ^ above are computed at the 
predictor step which has an error of the order of (^i)^, i e. the error (At)^ intro- 
duced in Eq. (6) is now eliminated. Therefore, with the expressions in Eqs. (15) and 
(16) one can use integration steps significantly larger than with the original Cundall 
spring. 

When considering discs, one does not take into account the shape of the par- 
ticles. Next, we consider the more realistic situation of irregular polygonal-shaped 
particles. Motion of rigid particles with polygonal shape is more complicated than 
that of simple discs, since the contact point no longer lies on the vector connecting 
the centers of mass. Further, for polygons, one must also be careful when decompos- 
ing the dynamics of each particle into translation and rotation around its center of 
mass. This implies recalculating each time the position of the center of mass (only 
from translation) and the relative position of the vertices (only from rotation). 

Therefore, for the translational contribution in Eqs. (15), we compute the 
overlap area between the two particles at t and at the predictor step. The overlap area 
is in general a polygon with a geometrical center that can be computed also at time t 
and the subsequent predictor step, yielding rdt) and respectively. The increment 
for the translational contribution will be just the tangential projection {tc (t) — r^) •t'^. 
Similarly, the contribution from the (polygonal) particle rotation is computed by 
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Fig. 10. Stress control test between two polygonal particles, as illustrated in Fig. 1. Com- 
parison of the relaxation time tn (in units of tc) when using the standard integration scheme 
(squares) and the proposed improved scheme (circles). Here, rotation is neglected. 



determining branch vectors, ri,{t) and r^, defined as the vectors joining the center 
of the particle and the center of the overlap area at time step t and the predictor 
step respectively. Computing the branch vector at t and at the predictor step, one 
derives the angle defined by them, namely 6 = arccos (r]^ • rf,(i)/(r^rf,(t))) and 
the average value (r^ — rf,(t))/2, yielding an increment in Eq. (16) given by 9{rf^ — 

n{t))/2. 

Figure 10 compares how the relaxation time varies with the normalized time 
step when the original Cundall approach is used (squares) and when our improved 
approach is introduced (circles). Clearly, the dependence on the integration step 
observed for the usual integration scheme disappears when our improved approach 
is introduced. Therefore, all the conclusions taken above for discs remain valid for 
polygons. 



5 Discussion and conclusions 

In this paper we introduced a technique to improve the accuracy of the numerical 
scheme used to compute the evolution of particle systems. 

To that end, we have first shown that the range of admissible integration steps has 
an upper limit significantly smaller than typically used. The accuracy of the numer- 
ical scheme not only depends on the associated error when computing the particle 
positions (predictor-corrector scheme), but also on the accuracy when determining 
the frictional force, which is usually implemented by the Cundall spring. Since the 
Cundall spring is linear in the integration step, the overall accuracy of the numerical 
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scheme cannot be higher than (At)'^. Therefore, when large integration steps are 
required, e.g. in slow shearing, the numerical scheme does not give accurate results. 

To overcome this problem we introduced an alternative approach for computing 
the frictional forces that suits not only the simple situation of discs but the more 
realistic situation of polygonal particles. Our approach is particularly suited for sit- 
uations where non-sliding contacts are relevant to the overall response. In general, 
for any other integration scheme, the substitution of the Cundall spring expression 
by the relations introduced in Eqs. (15) and (16), yields an error that is of the same 
order of the one associated with the predictor-corrector scheme. 

Inspired by the above results, some questions arise to further improve our ap- 
proach. First, the influence of the relations kt / kn and et /cn should also be consid- 
ered. Preliminary simulations have shown that the upper limit for the integration 
step increases with the value of kf/kn- Second, the test assumes a unique choice 
for the position of the contact point. However, in a system under shearing the inte- 
gration must be also performed before the appearance of new contacts. The initial 
contact point of a new contact will depend on the size of the integration step. This 
point should also be taken into account within our new approach, either by assuming 
some sort of interpolation or by using an event-driven scheme till the first contact 
point. Third, there is the problem of how to better define the contact point between 
two polygons. Since the contact point is taken as the geometrical center of their 
overlap area, the branch vectors also vary during rotation, which is not taken into 
account in our present approach. These and other points will be addressed in the 
future. 
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